Chapter 12 Source of functional differences

# A tibble: 4 × 2
  type_river count
  <chr>      <int>
1 absent         3
2 core         389
3 endemic       99
4 marginal      48

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

12.1 Genome size

# A tibble: 4 × 2
  type     genome_size
  <chr>          <dbl>
1 core        2835369.
2 endemic     2407507.
3 marginal    2396801.
4 <NA>        2030171 

# A tibble: 3 × 9
  .y.    group1  group2      n1    n2 statistic          p    p.adj p.adj.signif
* <chr>  <chr>   <chr>    <int> <int>     <dbl>      <dbl>    <dbl> <chr>       
1 length core    endemic    389    85     21740 0.00000533 0.000016 ****        
2 length core    marginal   389    62     15766 0.000101   0.000152 ***         
3 length endemic marginal    85    62      2738 0.688      0.688    ns          

12.2 Taxonomic variation

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  group_by(type,sample) %>% 
  summarise(fraction=sum(count)) %>% 
  group_by(type) %>% 
  summarise(mean(fraction))

`summarise()` has grouped output by 'type'. You can override using the `.groups` argument.

# A tibble: 4 × 2
  type     `mean(fraction)`
  <chr>               <dbl>
1 core               0.870 
2 endemic            0.0882
3 marginal           0.0407
4 <NA>               0.0746

12.3 Functional differences between fractions

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>%
    ggplot(aes(x=type, y=value, group=type))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type, p.adjust.method = "BH")

# A tibble: 3 × 10
  trait .y.   group1  group2      n1    n2 statistic          p      p.adj p.adj.signif
* <chr> <chr> <chr>   <chr>    <int> <int>     <dbl>      <dbl>      <dbl> <chr>       
1 mci   value core    endemic    389    85     21482 0.0000152  0.0000228  ****        
2 mci   value core    marginal   389    62     16678 0.00000126 0.00000378 ****        
3 mci   value endemic marginal    85    62      2903 0.294      0.294      ns          

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type, p.adjust.method = "BH")

# A tibble: 6 × 10
  trait        .y.   group1  group2      n1    n2 statistic           p       p.adj p.adj.signif
* <chr>        <chr> <chr>   <chr>    <int> <int>     <dbl>       <dbl>       <dbl> <chr>       
1 Biosynthesis value core    endemic    389    85    21545  0.0000118   0.0000177   ****        
2 Biosynthesis value core    marginal   389    62    17034. 0.000000179 0.000000537 ****        
3 Biosynthesis value endemic marginal    85    62     2986  0.169       0.169       ns          
4 Degradation  value core    endemic    389    85    20437  0.000644    0.002       **          
5 Degradation  value core    marginal   389    62    15153  0.001       0.002       **          
6 Degradation  value endemic marginal    85    62     2687  0.84        0.84        ns          

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>% 
    ggplot(aes(x=type, y=value, group=type))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type, p.adjust.method = "BH") %>% 
    print(n=100)

# A tibble: 42 × 10
   trait .y.   group1  group2      n1    n2 statistic           p       p.adj p.adj.signif
 * <chr> <chr> <chr>   <chr>    <int> <int>     <dbl>       <dbl>       <dbl> <chr>       
 1 B01   value core    endemic    389    85    20668. 0.000301    0.000451    ***         
 2 B01   value core    marginal   389    62    16226. 0.0000123   0.0000369   ****        
 3 B01   value endemic marginal    85    62     2914. 0.274       0.274       ns          
 4 B02   value core    endemic    389    85    21470. 0.0000159   0.0000238   ****        
 5 B02   value core    marginal   389    62    16529  0.00000274  0.00000822  ****        
 6 B02   value endemic marginal    85    62     2752. 0.649       0.649       ns          
 7 B03   value core    endemic    389    85    18336. 0.112       0.168       ns          
 8 B03   value core    marginal   389    62    13706. 0.081       0.168       ns          
 9 B03   value endemic marginal    85    62     2744. 0.666       0.666       ns          
10 B04   value core    endemic    389    85    16912. 0.74        0.74        ns          
11 B04   value core    marginal   389    62    15036  0.002       0.005       **          
12 B04   value endemic marginal    85    62     3186. 0.03        0.046       *           
13 B06   value core    endemic    389    85    22249  0.000000582 0.000000873 ****        
14 B06   value core    marginal   389    62    17100  0.000000123 0.000000369 ****        
15 B06   value endemic marginal    85    62     2915  0.273       0.273       ns          
16 B07   value core    endemic    389    85    22042  0.00000147  0.00000441  ****        
17 B07   value core    marginal   389    62    16465  0.0000038   0.0000057   ****        
18 B07   value endemic marginal    85    62     2798. 0.525       0.525       ns          
19 B08   value core    endemic    389    85    20472. 0.000563    0.000844    ***         
20 B08   value core    marginal   389    62    16182. 0.0000148   0.0000444   ****        
21 B08   value endemic marginal    85    62     3007  0.145       0.145       ns          
22 D01   value core    endemic    389    85    19124. 0.016       0.047       *           
23 D01   value core    marginal   389    62    12628. 0.527       0.527       ns          
24 D01   value endemic marginal    85    62     2368  0.238       0.357       ns          
25 D02   value core    endemic    389    85    21914. 0.00000254  0.00000762  ****        
26 D02   value core    marginal   389    62    15312. 0.000644    0.000966    ***         
27 D02   value endemic marginal    85    62     2430  0.422       0.422       ns          
28 D03   value core    endemic    389    85    19320. 0.015       0.022       *           
29 D03   value core    marginal   389    62    14923  0.003       0.008       **          
30 D03   value endemic marginal    85    62     2743  0.673       0.673       ns          
31 D05   value core    endemic    389    85    20434. 0.000649    0.000974    ***         
32 D05   value core    marginal   389    62    15667  0.000154    0.000462    ***         
33 D05   value endemic marginal    85    62     2868. 0.363       0.363       ns          
34 D06   value core    endemic    389    85    19654  0.006       0.018       *           
35 D06   value core    marginal   389    62    12642. 0.538       0.538       ns          
36 D06   value endemic marginal    85    62     2322  0.219       0.328       ns          
37 D07   value core    endemic    389    85    16382  0.895       0.972       ns          
38 D07   value core    marginal   389    62    12092. 0.972       0.972       ns          
39 D07   value endemic marginal    85    62     2690. 0.828       0.972       ns          
40 D09   value core    endemic    389    85    17899  0.225       0.225       ns          
41 D09   value core    marginal   389    62    15196. 0.000804    0.002       **          
42 D09   value endemic marginal    85    62     3028. 0.117       0.176       ns          

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>% 
    ggplot(aes(x=type, y=value, group=type))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

12.4 Functional differences between high and low endemisms

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment)) %>%
    ggplot(aes(x=type_prevalence_environment, y=value, group=type_prevalence_environment))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment)) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type_prevalence_environment, p.adjust.method = "BH")

# A tibble: 1 × 10
  trait .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif
* <chr> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 mci   value high   low       47    38       552 0.002 0.002 **          

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    ggplot(aes(x=type_prevalence_environment, y=value, group=type_prevalence_environment))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type_prevalence_environment, p.adjust.method = "BH") %>% 
    print(n=100)

# A tibble: 14 × 10
   trait .y.   group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
 * <chr> <chr> <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
 1 B01   value high   low       47    38      814. 0.485    0.485    ns          
 2 B02   value high   low       47    38      588. 0.007    0.007    **          
 3 B03   value high   low       47    38      668  0.045    0.045    *           
 4 B04   value high   low       47    38      684  0.065    0.065    ns          
 5 B06   value high   low       47    38      635  0.023    0.023    *           
 6 B07   value high   low       47    38      639  0.025    0.025    *           
 7 B08   value high   low       47    38      740. 0.179    0.179    ns          
 8 D01   value high   low       47    38      750. 0.143    0.143    ns          
 9 D02   value high   low       47    38      614  0.014    0.014    *           
10 D03   value high   low       47    38      592  0.008    0.008    **          
11 D05   value high   low       47    38      637  0.024    0.024    *           
12 D06   value high   low       47    38      682. 0.061    0.061    ns          
13 D07   value high   low       47    38      660. 0.04     0.04     *           
14 D09   value high   low       47    38      523  0.000929 0.000929 ***         

12.5 MAGs in different locations and shared among locations

12.5.1 Core

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="core") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

12.5.2 Endemic

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="endemic") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    ylim(0,1)+
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

12.5.3 Marginal

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="marginal") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    ylim(0,1)+
    scale_fill_manual(values=phylum_colors) +
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

12.5.4 Ordination

gift_pcoa <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame() %>%
    vegdist(method="euclidean") %>%
    pcoa()

gift_pcoa_rel_eigen <- gift_pcoa$values$Relative_eig[1:10]


# Get genome positions
gift_pcoa_vectors <- gift_pcoa$vectors %>% #extract vectors
  as.data.frame() %>% 
  select(Axis.1,Axis.2) # keep the first 2 axes

gift_pcoa_eigenvalues <- gift_pcoa$values$Eigenvalues[c(1,2)]

gift_pcoa_gifts <- cov(genome_gifts, scale(gift_pcoa_vectors)) %*% diag((gift_pcoa_eigenvalues/(nrow(genome_gifts)-1))^(-0.5)) %>%
  as.data.frame() %>% 
  rename(Axis.1=1,Axis.2=2) %>% 
  rownames_to_column(var="label") %>% 
  #get function summary vectors
  mutate(func=substr(label,1,3)) %>% 
  group_by(func) %>% 
  summarise(Axis.1=mean(Axis.1),
            Axis.2=mean(Axis.2)) %>% 
  rename(label=func) %>% 
  filter(!label %in% c("S01","S02","S03"))

set.seed(101)
scale <- 20 # scale for vector loadings
gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(genome_metadata, by="genome") %>%
  ggplot() +
      #genome positions
      scale_color_manual(values=phylum_colors)+
      geom_point(aes(x=Axis.1,y=Axis.2, color=phylum, size=length), 
                 alpha=0.9, shape=16) +
      #scale_color_manual(values=phylum_colors) +
      scale_size_continuous(range = c(0.1,5)) +
      #loading positions
      geom_segment(data=gift_pcoa_gifts, 
                   aes(x=0, y=0, xend=Axis.1 * scale, yend=Axis.2 * scale),
                    arrow = arrow(length = unit(0.3, "cm"), 
                    type = "open", 
                    angle = 25),
                    linewidth = 0.5, 
                    color = "black") +
     #Primary and secondary scale adjustments
     scale_x_continuous(name = paste0("PCoA1 (",round(gift_pcoa_rel_eigen[1]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA1")
            ) +
     scale_y_continuous(name = paste0("PCoA2 (",round(gift_pcoa_rel_eigen[2]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA2")
            ) +
    geom_label_repel(data = gift_pcoa_gifts,
                     aes(label = label, x = Axis.1 * scale, y = Axis.2 * scale),
                     segment.color = 'transparent') +
    xlim(-4.5,4.5) + 
    ylim(-3,2.5) +
    theme_minimal() + 
    theme(legend.position = "none")

structural <- core_microbiota %>% 
  mutate(type = case_when(
      (type_prevalence == "endemic") & (type_prevalence_environment == "high") ~ "present_high",
       (type_prevalence == "endemic") & (type_prevalence_environment == "low") ~ "present_low")) %>% 
    filter(!is.na(type)) %>% 
    select(genome,type)

enriched <- ancom_mag_gut_table %>% 
  mutate(type = case_when(
      lfc_environmentlow < 0 ~ "enriched_low",
      lfc_environmentlow >= 0 ~ "enriched_high")) %>% 
  select(genome,type)


scale <- 15 # scale for vector loadings
gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(bind_rows(structural,enriched), by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  mutate(type=ifelse(is.na(type),"neutral",type)) %>%
  group_by(type) %>%
  mutate(x_cen = median(Axis.1, na.rm = TRUE)) %>%
  mutate(y_cen = median(Axis.2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot() +
      scale_size_continuous(range = c(0.1,5)) +
      geom_point(aes(x=Axis.1,y=Axis.2, color=type, size=length), alpha=0.9, shape=16) +
      geom_segment(aes(x = x_cen, y = y_cen, xend = Axis.1, yend = Axis.2, color=type), alpha = 0.5) +
      scale_color_manual(values=c("#5A99D2","#F16144","#B7B7B7","#225072","#8C2C1C"))+
    theme_minimal() + 
    theme(legend.position = "none")

scale <- 15 # scale for vector loadings
gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(core_microbiota, by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  filter(type != "absent") %>% 
  mutate(shape = case_when(
      type == "core" ~ "core",
      type == "marginal" ~ "marginal",
      type == "endemic" & type_prevalence_environment == "high" ~ "endemic_high",
      type == "endemic" & type_prevalence_environment == "low" ~ "endemic_low"
    )) %>% 
  group_by(shape) %>%
  mutate(x_cen = median(Axis.1, na.rm = TRUE)) %>%
  mutate(y_cen = median(Axis.2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot() +
      #genome positions
      geom_segment(aes(x = x_cen, y = y_cen, xend = Axis.1, yend = Axis.2, color=shape), alpha = 0.5, show.legend = FALSE) +
      scale_color_manual(values=c("#9E8F71","#5A99D2","#F16144","#B7B7B7"))+
      theme_minimal() 

   # theme(legend.position = "none")

12.6 Functional variation

12.6.1 Core

# A tibble: 0 × 5
# ℹ 5 variables: gift <chr>, high <dbl>, low <dbl>, p_value <dbl>, p_adjust <dbl>

    Wilcoxon rank sum exact test

data:  mci by environment
W = 55, p-value = 0.01643
alternative hypothesis: true location shift is not equal to 0

# A tibble: 0 × 7
# ℹ 7 variables: gift <chr>, high <dbl>, low <dbl>, p_value <dbl>, p_adjust <dbl>, difference <dbl>, significance <chr>

term	df	SumOfSqs	R2	statistic	p.value
environment	1	0.7158447	0.06965384	2.136654	0.051
river	2	0.8505330	0.08275942	1.269336	0.223
Residual	26	8.7107971	0.84758674	NA	NA
Total	29	10.2771749	1.00000000	NA	NA

12.6.2 Endemic

    Wilcoxon rank sum exact test

data:  mci by environment
W = 24, p-value = 8.979e-05
alternative hypothesis: true location shift is not equal to 0

gift	high	low	p_value	p_adjust
B0104	0.318371850	0.477248051	4.266185e-03	8.738798e-03
B0106	0.652417204	0.750373702	2.349615e-02	3.825655e-02
B0204	0.272751511	0.411834564	9.967024e-04	2.751765e-03
B0205	0.349827863	0.607824698	2.150627e-04	7.803705e-04
B0207	0.349413509	0.607091324	8.978999e-05	4.072617e-04
B0208	0.277418167	0.571441654	1.125340e-04	4.610265e-04
B0209	0.418950756	0.654470027	1.665215e-03	3.916339e-03
B0210	0.314626667	0.599130873	2.701874e-03	5.815898e-03
B0211	0.505796266	0.763739750	2.514223e-06	6.450676e-05
B0212	0.341025004	0.546199073	4.939623e-03	9.802064e-03
B0213	0.415890582	0.550634314	1.643172e-02	2.782438e-02
B0215	0.255319614	0.554385824	5.777813e-04	1.834456e-03
B0216	0.105448827	0.426833768	2.041033e-05	1.993932e-04
B0217	0.281630700	0.400857555	6.955501e-04	2.154509e-03
B0220	0.108696049	0.075493630	9.874823e-03	1.717949e-02
B0221	0.255637379	0.383346961	9.874823e-03	1.717949e-02
B0303	0.168217937	0.258794974	2.643028e-04	8.833279e-04
B0309	0.028042752	0.096603777	7.399332e-03	1.408868e-02
B0310	0.000000000	0.039314500	8.257612e-03	1.519879e-02
B0402	0.556281087	0.493520756	1.408119e-03	3.804919e-03
B0601	0.382208213	0.577022449	2.654761e-05	2.408248e-04
B0602	0.506114317	0.745114749	1.403581e-04	5.401660e-04
B0603	0.260255972	0.467609340	5.618966e-05	3.161839e-04
B0701	0.432395583	0.679531792	8.978999e-05	4.072617e-04
B0703	0.022340365	0.156111414	5.726164e-05	3.161839e-04
B0704	0.312937738	0.598325101	3.432236e-05	2.564082e-04
B0705	0.184992803	0.398754859	1.554950e-05	1.645656e-04
B0706	0.336719078	0.509754909	7.543545e-03	1.408868e-02
B0707	0.584533168	0.687871927	7.543545e-03	1.408868e-02
B0710	0.016909779	0.108809582	1.581108e-03	3.916339e-03
B0711	0.208154526	0.416638817	7.121052e-05	3.617495e-04
B0712	0.099046619	0.204586247	5.703508e-03	1.114378e-02
B0804	0.456763954	0.763616891	1.177172e-05	1.359099e-04
B0805	0.042739966	0.173618219	5.802053e-07	4.360009e-05
B0901	0.100157586	0.041913672	1.665215e-03	3.916339e-03
B0903	0.000000000	0.019080084	2.539636e-06	6.450676e-05
B1004	0.125149306	0.140982835	8.964396e-03	1.603490e-02
B1012	0.001356128	0.019240248	1.619167e-04	6.048064e-04
B1014	0.022601362	0.000000000	2.539636e-06	6.450676e-05
D0104	0.058683884	0.159932473	2.306458e-03	5.050347e-03
D0201	0.092234219	0.260028653	8.340902e-04	2.407488e-03
D0202	0.160795013	0.371123533	5.618966e-05	3.161839e-04
D0203	0.342605111	0.471551663	3.432236e-05	2.564082e-04
D0204	0.250243216	0.395633777	2.306458e-03	5.050347e-03
D0205	0.063752861	0.169920437	7.121052e-05	3.617495e-04
D0206	0.147551315	0.443090128	5.618966e-05	3.161839e-04
D0207	0.412156668	0.558207530	2.090231e-02	3.447523e-02
D0208	0.155823369	0.301738922	5.618966e-05	3.161839e-04
D0209	0.131997770	0.352678185	6.549873e-06	9.242598e-05
D0210	0.133702102	0.284613983	1.177172e-05	1.359099e-04
D0212	0.133552981	0.422562064	6.549873e-06	9.242598e-05
D0213	0.126218478	0.306630447	8.340902e-04	2.407488e-03
D0305	0.366492686	0.423996527	4.939623e-03	9.802064e-03
D0306	0.146105878	0.434641186	6.549873e-06	9.242598e-05
D0501	0.790847790	0.872992541	1.805380e-02	3.016885e-02
D0502	0.480358478	0.285918054	4.781407e-04	1.557022e-03
D0503	0.125879327	0.000000000	6.866156e-07	4.360009e-05
D0505	0.286597240	0.350858031	8.641809e-03	1.567871e-02
D0506	0.211360464	0.391670686	6.549873e-06	9.242598e-05
D0507	0.110052668	0.210004157	1.125340e-04	4.610265e-04
D0509	0.440041988	0.526487614	1.125601e-02	1.931774e-02
D0510	0.000000000	0.013573683	2.308808e-04	7.924829e-04
D0511	0.611751724	0.767405683	1.665215e-03	3.916339e-03
D0513	0.224858211	0.438902945	3.155375e-03	6.678878e-03
D0601	0.034174614	0.141369358	2.997491e-02	4.758517e-02
D0602	0.000000000	0.030499587	2.308808e-04	7.924829e-04
D0606	0.003034342	0.043637689	5.152105e-05	3.161839e-04
D0607	0.008552574	0.041697052	2.096373e-03	4.840715e-03
D0608	0.000000000	0.002444638	3.625355e-03	7.547870e-03
D0609	0.158687978	0.084565489	1.403581e-04	5.401660e-04
D0610	0.000000000	0.007160308	1.522953e-03	3.916339e-03
D0702	0.455688142	0.326464654	2.635404e-02	4.236663e-02
D0705	0.323241782	0.258159545	2.306458e-03	5.050347e-03
D0706	0.059976159	0.148351231	1.665215e-03	3.916339e-03
D0807	0.015074471	0.041747020	8.978999e-05	4.072617e-04
D0816	0.043641253	0.091390659	8.340902e-04	2.407488e-03
D0817	0.005827776	0.019925548	1.445050e-03	3.823361e-03
D0907	0.162914556	0.328681837	9.967024e-04	2.751765e-03
D0908	0.305403136	0.623522292	3.432236e-05	2.564082e-04
D0910	0.124407773	0.396891803	1.125340e-04	4.610265e-04

gift	high	low	p_value	p_adjust	difference	significance
B0104	0.318371850	0.477248051	4.266185e-03	8.738798e-03	-0.158876201	significant
B0106	0.652417204	0.750373702	2.349615e-02	3.825655e-02	-0.097956497	significant
B0204	0.272751511	0.411834564	9.967024e-04	2.751765e-03	-0.139083053	significant
B0205	0.349827863	0.607824698	2.150627e-04	7.803705e-04	-0.257996835	significant
B0207	0.349413509	0.607091324	8.978999e-05	4.072617e-04	-0.257677815	significant
B0208	0.277418167	0.571441654	1.125340e-04	4.610265e-04	-0.294023487	significant
B0209	0.418950756	0.654470027	1.665215e-03	3.916339e-03	-0.235519272	significant
B0210	0.314626667	0.599130873	2.701874e-03	5.815898e-03	-0.284504206	significant
B0211	0.505796266	0.763739750	2.514223e-06	6.450676e-05	-0.257943484	significant
B0212	0.341025004	0.546199073	4.939623e-03	9.802064e-03	-0.205174069	significant
B0213	0.415890582	0.550634314	1.643172e-02	2.782438e-02	-0.134743732	significant
B0215	0.255319614	0.554385824	5.777813e-04	1.834456e-03	-0.299066210	significant
B0216	0.105448827	0.426833768	2.041033e-05	1.993932e-04	-0.321384940	significant
B0217	0.281630700	0.400857555	6.955501e-04	2.154509e-03	-0.119226854	significant
B0220	0.108696049	0.075493630	9.874823e-03	1.717949e-02	0.033202420	significant
B0221	0.255637379	0.383346961	9.874823e-03	1.717949e-02	-0.127709582	significant
B0303	0.168217937	0.258794974	2.643028e-04	8.833279e-04	-0.090577037	significant
B0309	0.028042752	0.096603777	7.399332e-03	1.408868e-02	-0.068561026	significant
B0310	0.000000000	0.039314500	8.257612e-03	1.519879e-02	-0.039314500	significant
B0402	0.556281087	0.493520756	1.408119e-03	3.804919e-03	0.062760331	significant
B0601	0.382208213	0.577022449	2.654761e-05	2.408248e-04	-0.194814235	significant
B0602	0.506114317	0.745114749	1.403581e-04	5.401660e-04	-0.239000432	significant
B0603	0.260255972	0.467609340	5.618966e-05	3.161839e-04	-0.207353368	significant
B0701	0.432395583	0.679531792	8.978999e-05	4.072617e-04	-0.247136209	significant
B0703	0.022340365	0.156111414	5.726164e-05	3.161839e-04	-0.133771050	significant
B0704	0.312937738	0.598325101	3.432236e-05	2.564082e-04	-0.285387364	significant
B0705	0.184992803	0.398754859	1.554950e-05	1.645656e-04	-0.213762056	significant
B0706	0.336719078	0.509754909	7.543545e-03	1.408868e-02	-0.173035831	significant
B0707	0.584533168	0.687871927	7.543545e-03	1.408868e-02	-0.103338760	significant
B0710	0.016909779	0.108809582	1.581108e-03	3.916339e-03	-0.091899802	significant
B0711	0.208154526	0.416638817	7.121052e-05	3.617495e-04	-0.208484291	significant
B0712	0.099046619	0.204586247	5.703508e-03	1.114378e-02	-0.105539628	significant
B0804	0.456763954	0.763616891	1.177172e-05	1.359099e-04	-0.306852936	significant
B0805	0.042739966	0.173618219	5.802053e-07	4.360009e-05	-0.130878253	significant
B0901	0.100157586	0.041913672	1.665215e-03	3.916339e-03	0.058243913	significant
B0903	0.000000000	0.019080084	2.539636e-06	6.450676e-05	-0.019080084	significant
B1004	0.125149306	0.140982835	8.964396e-03	1.603490e-02	-0.015833529	significant
B1012	0.001356128	0.019240248	1.619167e-04	6.048064e-04	-0.017884120	significant
B1014	0.022601362	0.000000000	2.539636e-06	6.450676e-05	0.022601362	significant
D0104	0.058683884	0.159932473	2.306458e-03	5.050347e-03	-0.101248589	significant
D0201	0.092234219	0.260028653	8.340902e-04	2.407488e-03	-0.167794434	significant
D0202	0.160795013	0.371123533	5.618966e-05	3.161839e-04	-0.210328521	significant
D0203	0.342605111	0.471551663	3.432236e-05	2.564082e-04	-0.128946551	significant
D0204	0.250243216	0.395633777	2.306458e-03	5.050347e-03	-0.145390561	significant
D0205	0.063752861	0.169920437	7.121052e-05	3.617495e-04	-0.106167576	significant
D0206	0.147551315	0.443090128	5.618966e-05	3.161839e-04	-0.295538813	significant
D0207	0.412156668	0.558207530	2.090231e-02	3.447523e-02	-0.146050862	significant
D0208	0.155823369	0.301738922	5.618966e-05	3.161839e-04	-0.145915553	significant
D0209	0.131997770	0.352678185	6.549873e-06	9.242598e-05	-0.220680415	significant
D0210	0.133702102	0.284613983	1.177172e-05	1.359099e-04	-0.150911881	significant
D0212	0.133552981	0.422562064	6.549873e-06	9.242598e-05	-0.289009083	significant
D0213	0.126218478	0.306630447	8.340902e-04	2.407488e-03	-0.180411970	significant
D0305	0.366492686	0.423996527	4.939623e-03	9.802064e-03	-0.057503841	significant
D0306	0.146105878	0.434641186	6.549873e-06	9.242598e-05	-0.288535308	significant
D0501	0.790847790	0.872992541	1.805380e-02	3.016885e-02	-0.082144750	significant
D0502	0.480358478	0.285918054	4.781407e-04	1.557022e-03	0.194440424	significant
D0503	0.125879327	0.000000000	6.866156e-07	4.360009e-05	0.125879327	significant
D0505	0.286597240	0.350858031	8.641809e-03	1.567871e-02	-0.064260791	significant
D0506	0.211360464	0.391670686	6.549873e-06	9.242598e-05	-0.180310222	significant
D0507	0.110052668	0.210004157	1.125340e-04	4.610265e-04	-0.099951489	significant
D0509	0.440041988	0.526487614	1.125601e-02	1.931774e-02	-0.086445625	significant
D0510	0.000000000	0.013573683	2.308808e-04	7.924829e-04	-0.013573683	significant
D0511	0.611751724	0.767405683	1.665215e-03	3.916339e-03	-0.155653959	significant
D0513	0.224858211	0.438902945	3.155375e-03	6.678878e-03	-0.214044734	significant
D0601	0.034174614	0.141369358	2.997491e-02	4.758517e-02	-0.107194744	significant
D0602	0.000000000	0.030499587	2.308808e-04	7.924829e-04	-0.030499587	significant
D0606	0.003034342	0.043637689	5.152105e-05	3.161839e-04	-0.040603346	significant
D0607	0.008552574	0.041697052	2.096373e-03	4.840715e-03	-0.033144478	significant
D0608	0.000000000	0.002444638	3.625355e-03	7.547870e-03	-0.002444638	significant
D0609	0.158687978	0.084565489	1.403581e-04	5.401660e-04	0.074122489	significant
D0610	0.000000000	0.007160308	1.522953e-03	3.916339e-03	-0.007160308	significant
D0702	0.455688142	0.326464654	2.635404e-02	4.236663e-02	0.129223488	significant
D0705	0.323241782	0.258159545	2.306458e-03	5.050347e-03	0.065082237	significant
D0706	0.059976159	0.148351231	1.665215e-03	3.916339e-03	-0.088375073	significant
D0807	0.015074471	0.041747020	8.978999e-05	4.072617e-04	-0.026672549	significant
D0816	0.043641253	0.091390659	8.340902e-04	2.407488e-03	-0.047749406	significant
D0817	0.005827776	0.019925548	1.445050e-03	3.823361e-03	-0.014097772	significant
D0907	0.162914556	0.328681837	9.967024e-04	2.751765e-03	-0.165767281	significant
D0908	0.305403136	0.623522292	3.432236e-05	2.564082e-04	-0.318119156	significant
D0910	0.124407773	0.396891803	1.125340e-04	4.610265e-04	-0.272484030	significant

term	df	SumOfSqs	R2	statistic	p.value
environment	1	18.732917	0.27521376	11.099093	0.001
river	2	5.451379	0.08008868	1.614948	0.148
Residual	26	43.882492	0.64469756	NA	NA
Total	29	68.066788	1.00000000	NA	NA

12.6.3 Marginal

    Wilcoxon rank sum exact test

data:  mci by environment
W = 56, p-value = 0.03277
alternative hypothesis: true location shift is not equal to 0

# A tibble: 33 × 5
   gift     high    low  p_value p_adjust
   <chr>   <dbl>  <dbl>    <dbl>    <dbl>
 1 B0221 0.177   0.527  0.00268    0.0297
 2 B0605 0.0886  0.396  0.00663    0.0371
 3 B0709 0       0.0329 0.000980   0.0196
 4 B1028 0.00776 0.0891 0.00224    0.0297
 5 D0101 0.00637 0.0634 0.00261    0.0297
 6 D0103 0.0689  0.216  0.00839    0.0435
 7 D0201 0.0589  0.252  0.00276    0.0297
 8 D0202 0.110   0.323  0.00820    0.0435
 9 D0203 0.214   0.424  0.00932    0.0450
10 D0204 0.225   0.483  0.00631    0.0368
# ℹ 23 more rows

gift	high	low	p_value	p_adjust	difference	significance
B0221	0.177144655	0.52701883	2.677235e-03	0.029723307	-0.34987417	significant
B0605	0.088604745	0.39623709	6.627165e-03	0.037112124	-0.30763234	significant
B0709	0.000000000	0.03290746	9.803796e-04	0.019607591	-0.03290746	significant
B1028	0.007758061	0.08906279	2.236734e-03	0.029723307	-0.08130473	significant
D0101	0.006371054	0.06336166	2.607633e-03	0.029723307	-0.05699060	significant
D0103	0.068861705	0.21552103	8.387951e-03	0.043493078	-0.14665933	significant
D0201	0.058876923	0.25212032	2.760021e-03	0.029723307	-0.19324339	significant
D0202	0.110126355	0.32295102	8.200669e-03	0.043493078	-0.21282467	significant
D0203	0.213841081	0.42380897	9.322462e-03	0.045004989	-0.20996789	significant
D0204	0.225490303	0.48319440	6.306083e-03	0.036785487	-0.25770410	significant
D0205	0.064765347	0.17098720	9.226787e-03	0.045004989	-0.10622185	significant
D0208	0.113400443	0.25966177	4.235560e-03	0.032943248	-0.14626132	significant
D0209	0.101669651	0.28432354	1.758623e-03	0.027356358	-0.18265389	significant
D0212	0.172471768	0.29920957	1.048346e-02	0.046042591	-0.12673780	significant
D0213	0.094273805	0.26659680	3.217419e-03	0.030029244	-0.17232300	significant
D0301	0.000000000	0.07340385	6.096424e-03	0.036785487	-0.07340385	significant
D0302	0.004398428	0.08700617	2.275440e-04	0.006371233	-0.08260774	significant
D0304	0.009925215	0.22320473	5.606955e-05	0.002616579	-0.21327951	significant
D0307	0.015936692	0.24428734	1.718478e-04	0.006014672	-0.22835064	significant
D0308	0.119232778	0.32375506	5.354418e-03	0.036785487	-0.20452228	significant
D0502	0.059489928	0.55204220	1.708261e-05	0.001195782	-0.49255228	significant
D0503	0.018551988	0.14195916	1.085290e-02	0.046042591	-0.12340717	significant
D0509	0.273367616	0.53084131	1.054085e-02	0.046042591	-0.25747370	significant
D0518	0.092542585	0.27773775	4.091025e-03	0.032943248	-0.18519517	significant
D0603	0.000000000	0.10245182	6.096424e-03	0.036785487	-0.10245182	significant
D0607	0.009275994	0.10179302	6.834773e-04	0.015947804	-0.09251702	significant
D0611	0.000000000	0.14970317	6.096424e-03	0.036785487	-0.14970317	significant
D0704	0.163489399	0.43901684	3.994756e-03	0.032943248	-0.27552744	significant
D0706	0.003315462	0.10115430	5.202430e-03	0.036785487	-0.09783884	significant
D0807	0.006303428	0.21451084	1.251696e-05	0.001195782	-0.20820742	significant
D0816	0.035535517	0.22732398	3.217419e-03	0.030029244	-0.19178846	significant
D0817	0.002166158	0.01700964	1.492190e-03	0.026113324	-0.01484348	significant
D0904	0.095233051	0.00000000	1.053795e-02	0.046042591	0.09523305	significant

term	df	SumOfSqs	R2	statistic	p.value
environment	1	22.66237	0.11067970	3.415318	0.021
river	2	16.20626	0.07914899	1.221177	0.235
Residual	25	165.88774	0.81017131	NA	NA
Total	28	204.75637	1.00000000	NA	NA